##

rm(list = ls())

## 

library(tidyverse)
library(lfe)
library(lubridate)
library(broom)

source('code/function_tidy_felm.R')

## Get data
## Note that this table requires data that we cannot make available publicly

civ_mip <- read_rds("data_private/civey_issue_replication.rds")

## Monthly intervals

start <- as.Date("2021-01-01")
end <- as.Date("2021-10-01")

dates_m <- floor_date(seq(start, end, by = 'month'), 
                      unit = "month") + 11

## Add time periods

civ_mip <- civ_mip %>% 
  mutate(month_10 = cut(date, dates_m)) %>% 
  mutate(month_10 = as.numeric(as.factor(month_10)) - 7) %>% 
  filter(!is.na(month_10)) %>% 
  mutate(post = ifelse(month_10 > 0, 1, 0))

## Select DVs

dv_list <- c("climate",
             "interior",
             "economy",
             "migr",
             "foreign_pol",
             "soc_pol")

## Loop to run separate models for each outcome

out <- dv_list %>% 
  lapply(function(dv) {
    
    f1 <- paste0(dv, " ~ treat_heavy*post + 
             treat_weakly*post |ags + month_10 + gender + educ + age|0|id + ags") %>% 
      as.formula()
    
    m <- felm(f1, data = civ_mip)
    
  })

## Mean of DV

dv_mean <- out %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == 'treat_heavy') %>% 
  pull(dv_mean)%>% 
  round(3)

## SD of DV

dv_sd <- out %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == 'treat_heavy') %>% 
  pull(dv_sd) %>% 
  round(3)

#### Table A.3 - effect on other issues ####

stargazer::stargazer(out, keep = 'post:treat_weakly|treat_heavy:post',
                     order = c("post:treat_weakly", "treat_heavy:post"),
                     covariate.labels = c('Weakly affected * post flood',
                                          'Highly affected * post flood'),
                     style = 'ajps',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_mean),
                                      c('DV SD', dv_sd),
                                      c("Month FE", rep("Yes", 6)),
                                      c("County FE", rep("Yes", 6))))
